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Abstract — The migration of an isolated gas bubble in an immiscible liquid possessing a temperature 
gradient is analyzed in the absence of gravity. The driving force for the bubble motion is the shear stress 
at the interface which is a consequence of the temperature dependence of the surface tension. The analysis 
is performed under conditions for which the Marangoni number is large, i.e. energy is transferred 
predominantly by convection. Velocity fields in the limit of both small and large Reynolds numbers are 
used. The thermal problem is treated by standard boundary layer theory. The outer temperature field is 
obtained in the vicinity of the bubble. A similarity solution is obtained for the inner temperature field. 
For both small and large Reynolds numbers, the asymptotic values of the scaled migration velocity of 
the bubble in the limit of large Marangoni numbers are calculated. The results show that the migration 
velocity has the same scaling for both low and large Reynolds numbers, but with a different coefficient. 
Higher order thermal boundary layers are analyzed for the large Reynolds number flow field and the 
higher order corrections to the migration velocity are obtained. Results are also presented for the 
momentum boundary layer and the thermal wake behind the bubble, for large Reynolds number 
conditions. Copyright © 1996 Elsevier Science Ltd. 

Key Words: thermocapillary migration, bubbles and drops, Marangoni convection, boundary layers, 
asymptotic analysis 


1. INTRODUCTION 

With the advent of the space program and the potential for space processing, there has been 
considerable interest in the motion of bubbles and drops in reduced gravity, both for fundamental 
understanding and actual applications (see review papers cited below). Thermally-induced surface 
tension gradient driven flow is the most often studied (and encountered) mechanism for this 
motion. When a bubble or drop is placed in a continuous phase in which a temperature gradient 
exists, the variation of interfacial tension with temperature along the interface will lead to a 
tangential stress. This stress will cause the motion of neighboring fluid and will lead to propulsion 
of the bubble in the direction of warmer fluid since surface tension decreases with temperature in 
most common fluids. Such motion is termed “thermocapillary migration” and forms the subject 
of the present study. 

The two most important parameters governing thermocapillary motion are the Reynolds and 
Marangoni numbers. The Reynolds number is R q=V r RJv and the Marangoni number is 
Ma = V R RJa. Here R l is the radius of the bubble, v and a are the kinematic viscosity and the 
thermal diffusivity of the surrounding liquid. The characteristic velocity V R = (-a r )ARJ pv where 
a denotes the surface tension, <r T is the rate of change of surface tension with temperature (assumed 
to be a constant, typically negative, in this analysis) and A is the temperature gradient in the liquid. 
In thermocapillary migration, the Marangoni number serves the role of a Peclet number, its 
magnitude signifying the relative importance of convective transport of energy when compared to 
conduction. 

Young et al. (1959) performed pioneering theoretical and experimental work on air bubbles in 
a silicone oil under conditions of negligible Reynolds and Marangoni numbers. Much has been 
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done subsequently and most of the existing literature has been reviewed in Wozniak et al. (1988) 
and Subramanian (1992). Bratukhin (1976) and Thompson et al. (1980) attempted to improve upon 
the solution of Young et al. by a regular perturbation scheme in the Reynolds number, so as to 
include effects of inertia and shape deformations. Subramanian (1981) pointed out that the above 
regular perturbation approach is unable to satisfy the far-field condition on the temperature field 
at higher orders and obtained a solution of the energy equation for small values of the Marangoni 
number by the method of matched asymptotic expansions. Crespo & Manuel (1983) and 
Balasubramaniam & Chai (1987) showed that the solution of Young et al. is an exact solution of 
the momentum equation for any Reynolds number, so long as the Marangoni number is negligible. 
The latter authors also calculated the small deformations of a drop from a spherical shape. 
Numerical calculations of the velocity and temperature fields for a spherical bubble as well as the 
bubble velocity have been performed by Szymczyk & Siekmann (1988), Shankar & Subramanian 
(1988) and Balasubramaniam & Lavery (1989), and the deformation of a bubble and its impact 
on the migration velocity have been calculated by Haj-Hariri et al. (1990), Chen & Lee (1992) and 
Nas & Tryggvason (1993). The transient migration of a drop from an initial state of rest to the 
steady state has been analyzed by Dill & Balasubramaniam (1992). 

The objective of the present study is to analyze the motion of an isolated gas bubble subjected 
to a uniform temperature gradient under zero gravity conditions when the Marangoni number is 
large. We consider two limiting cases, Re->0 and Re->co. Physically, the former is the case of a 
highly viscous flow of a large Prandtl number fluid (Pr = v/a). The latter corresponds to an 
inertia-dominated flow of a fluid of low to moderate Prandtl number. An example of the first is 
the movement of a bubble of radius 5 mm in a silicone oil (which is a model fluid used in many 
terrestrial and space-flight experiments) of nominal kinematic viscosity 50 centistokes that is 
subjected to a temperature gradient of 1 K/mm. The Reynolds number in this case is around 0.5 
and the Marangoni number is around 240. An example for the small Prandtl number case is the 
motion of a bubble of radius 2 mm in a melt of the semiconductor silicon with a temperature 
gradient of 2.5 K/mm. Here the Reynolds number is around 1700 and the Marangoni number is 
roughly 50. 

It is assumed that the bubble-liquid interface is clean and free from surfactants. Several such 
interfaces have been experimented with in the context of surface-driven motion and bubble 
migration, notably the air-silicone oil interface. Bubbles are also modeled as undeformed spheres 
in this analysis. In the large Pr example given above, the Capillary number (Ca = pvV R / a) is 0.01 
and deformation is negligible. In the small Pr example the Capillary number is 0.003, but the Weber 
number (ReCa) is not small and deformations may not be neglected. However experimental and 
numerical evidence thus far suggests that deformation is negligible and nowhere near what is 
observed for gravity-driven motion. Nas & Tryggvason (1993) performed a two-dimensional 
calculation with Re = 2000, Ma = 400 and Ca = 0.0166 and report negligible deformations. Also 
no flow separation is evident even at Re = 2000 in the calculations by Nas & Tryggvason and 
Balasubramaniam & Lavery. 

To keep the analysis tractable it is assumed that viscosity is not a function of temperature. While 
viscosity changes due to temperature variations can be important in virtually all experiments, a 
largely successful approach in prediction of the bubble velocity has been to consider such variation 
“quasi-static”, i.e. a local value of the viscosity is used in an expression obtained from a theory 
which assumes constant viscosity. 

In cases of both small and large Re, for large Ma, convective energy transport dominates over 
conduction almost everywhere. It will be seen that a thermal boundary layer must exist near the 
bubble surface in order to satisfy the boundary conditions. In the analogous fluid mechanical 
problem in the limit Re-»oo (Chao 1962; Moore 1963), the boundary layer results in 0(1) change 
in the stress; the velocity changes in the boundary layer are small, of 0(Re _1/2 ). In constrast, it 
will be shown that the thermal boundary layer produces temperature field changes of 0(1) here. 

A key question is the dependence of the scaled bubble velocity on Ma. From their numerical 
results for the temperature field under negligible Re, Shankar & Subramanian (1988) conclude that 
the scaled bubble velocity v x (non-dimensionalized by the velocity scale K R ) shows a logarithmi- 
cally decreasing behavior with Ma in the range 20 < Ma < 200. A fit to the numerical results of 
the form ^ = 1.59/(1.84 + In Ma) was found to be indistinguishable from the numerical results in 
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that range. The extrapolated value of from this fit goes to zero as Ma becomes infinite. 
Balasubramaniam & La very (1989) performed numerical computations in the range 
0 ^ Ma < 1000, 0 ^ Re ^ 2000. Their results show that while v m decreased with increasing Ma at 
fixed Re in this range, no clear asymptotic behavior could be discerned. 

While the present analysis was under way, Crespo & Jimenez-Fernandez (1991, 1992) have 
performed thermal boundary layer analyses for large Ma using velocity fields applicable in the 
limits of both small and large Re. Their analysis is very similar to ours presented in section 3. It 
will be shown that their result for the large Re case is in agreement with the present solution. 
However, their result for the Re->0 case appears to be in error in the numerical calculations for 
the Fourier coefficients in the series expansion of the velocity field. Also, in the case of large Re, 
we present the higher order corrections to the bubble migration velocity. 

In what follows, section 2 summarizes the known velocity fields for small and large Re, section 
3 deals with the boundary layer formalism for the leading order temperature field and section 4 
deals with the higher order thermal boundary layers for the large Re flow field. The results are 
discussed in section 5 and some concluding remarks are made in section 6. The momentum 
boundary layer for large Re is discussed in appendix A and the thermal wake for large Reynolds 
numbers is analyzed in appendix B. 


2. VELOCITY FIELDS FOR SMALL AND LARGE REYNOLDS NUMBERS 


Subramanian (1981) has analyzed the thermocapillary migration of a gas bubble in the limit of 
zero Reynolds number and small Marangoni numbers. The general solution for the flow field 
appropriate to axisymmetric thermocapillary migration problems in a reference frame traveling 
with the bubble is reproduced below from that analysis (see also Levan & Newman 1976). Thus 
the scaled flow variables in the limit of small Re are given by 


<Hr, = 7(1 — - rj + ^t n(n - I)/,/— - 

u(r, n) = \ /*(l - p) + j t »(« - 1 V.(^rr ■ - 

T ( 1 \ 1 00 (n — 3 

v(r, n) = —j (1 - fi 2 ) ll2 y 2 + + ~ 5>(« - + 


1 ) 

1 cm 

[i] 

r-V 

■W, 

(i u) 

[2] 

1 — n\ 

i cm 

[3] 

r n+\ ) 

'(\-n 2 ) 112 


Lengths are scaled by the bubble radius R x and velocities by V R defined earlier. The scaled radial 
coordinate is r and n = cos 6. The polar angle 6 is measured from the front stagnation streamline. 
The scaled radial and tangential velocities are u and v, respectively, and t ft is the scaled 
streamfunction. C n {fx) and P„(/x) are the Gegenbauer function of order n and degree —1/2 and 
the Legendre polynomial of order n, respectively. The Fourier coefficients /„ are obtained from the 
temperature field on the bubble surface: 


f 1 dT 

I n =~ C n (pi)—(l,n)dn n^2 

J-i 


dpi 


[4] 


Here, f is a transformed temperature related to the physical temperature t by 
T = {f — Av a, V R t)/(AR l ), where A is the temperature gradient far away from the bubble and t is 
time. The scaled bubble migration velocity is 


[5] 

which is obtained from the condition that the net hydrodynamic force on the migrating bubble is 
zero. 

When the Reynolds number is large, the velocity field may be approximated by a potential flow 
field. As in body force driven motion (Levich 1962), there is a thin boundary layer near the bubble 
surface in the stress field in order to satisfy the tangential stress balance at the surface (unlike in 
body force-driven motion, the surface shear stress in non-zero and equals the surface tension 
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gradient). However, the velocity field, to leading order, is that corresponding to potential flow all 
the way to the bubble surface. The correction to this field in the boundary layer occurs at 0(Re ' 1/2 ) 
as discussed in appendix A, and is neglected in the present analysis. 

The potential flow field does not exert any force on the bubble regardless of the bubble velocity. 
Thus an alternate condition is required to calculate v w . This is provided by viscous dissipation 
arguments (Levich 1962; Crespo & Manuel 1983; Balasubramaniam 1987)— at steady state 
migration, the rate at which work is done by the surface tension force equals the rate at which 
energy is dissipated by viscosity. 

Thus the large Reynolds number flow field, again in a reference frame traveling with the bubble, 
is given below 

ip(r,9) = y sin 2 0^r 2 -i^ u(r,9) = -v m cos 0^1 

^,/O = y«>sin0^1 + t 6 3 

The scaled bubble velocity is 

where the latter is obtained by the dissipation argument mentioned above. It may be noticed that 
[6] and [7] for the large Re case may be obtained from [1], [4] and [5] for the small Re case formally 
by setting I n - 0, n ^ 3. Thus the large Re problem may be treated as a special case of the small 
Re problem in this analysis. 


3. LEADING ORDER TEMPERATURE FIELD 

The temperature gradient far away from the bubble is assumed to be a constant. Typically, the 
temperature coefficient of surface tension, do jdT, is negative. The bubble would then move toward 
the warm portion of the surrounding liquid. In a reference frame moving with the bubble, the 
temperature field is unsteady because the field at infinity is unsteady in this reference frame. 
However, a modified temperature field obtained by subtracting from it the undisturbed temperature 
field evaluated at the center of the bubble will approach a steady field at infinity and therefore will 
achieve a steady state (Subramanian 1981). Thus the energy equation for this steady field T(r, 9) 


may be written as (Balasubramaniam & Chai 1987) 

8T , vdT 1 ^ 
00 + dr r 39 Ma 

[8] 

with boundary conditions 

-^(i./O-o 

[9] 

i.e. the bubble surface is adiabatic, 

T-+rfi as r->oo 

[10] 


i.e. the temperature field approaches the undisturbed field far away. In [8], u and v are given by 
[2] and [3] and is obtained from [5]. 

Equations [8]— [10], for known v nt are linear. Nevertheless analytical solution is formidable and 
we resort to perturbation theory for solution in the limit Ma-*co. For well-known reasons, the 
perturbation is singular, and therefore we use the method of matched asymptotic expansions. For 
large Ma, one would expect a radial thermal boundary layer to be present near the bubble surface. 
As usual, the temperature field outside the boundary layer is called the outer temperature field and 
that within the boundary layer is termed the inner temperature field. 
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3.1. Outer temperature field 
Let 



I. = I M + o( 1) n» 2 [12] 

u — u 0 + • • • [13] 

v=v 0 + [14] 
+ * * • [15] 

r(r,//) = T 0 (r,p) + • • • [16] 


where T is the outer temperature field. Substituting [13]-[16] in [8], the energy equation to leading 
order is 


dT 0 u 0 2 ^ 1/2 3T 0 . 

v x o + u 0 - (1 -^ 2 y i2 — = 0 

or r op, 

T 0 -*rp as r->oo 


[17] 


[18] 


Transforming from (r, p) to (r, t p 0 ) coordinates where ip is the streamfunction, the solution may 
be written as 


T 0 (r,p) = rp + 


v 1 

— ( — ^co 0 + y o s l n ^ ~ M o cos0)d ? 

, oo ^0 


[19] 


where the integrand must be expressed as a function of (r, i p 0 ). This may be further specialized for 
the large Re flow field to be 


T 0 (r, p) = rp + 


1 


3<Ao 
^ooo^ 2 - I/O 


- 1 


( 0 - 1 ) 


1 - 


2ip 0 


fcoo(0 - 1/0 


1/2 


dr 


[20] 


We are unable to find a closed form result for the integrals in [19] or [20]. Specialized results may 
be obtained for (i) r » 1 and (ii) 0 = 0, n for large Re. As can be expected, the outer solution is 
unable to satisfy the boundary conditions at the bubble surface and a thermal boundary layer exists 
there. To analyze this boundary layer and perform the necessary matching, the outer solution is 
needed near the bubble surface. We obtain it below via an expansion of the velocities near r = 1. 

3. 1. 1. Outer field near r — 1 . Near r = 1 , the leading order results for the velocities may be 
expanded as 


u a (r,n) = \l^(r — t »(" - l)I„,P,-M+0((r - l) 2 ) [21] 

Z Z n = 3 

3 1 00 C (u) 

v 0 (r, p) = -- (1 - p 2 y /2 I 2o - j S «(» - l) / «o' (if~^ + °( r ~ t 22 l 

Using these in [17], the solution for T 0 can be obtained as 

T 0 (r,p) = ^-\n(r - 1) -h In (f(p)) + g(ip) [23] 

A 0 


A 0 = 3 — 


[24] 
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[ 25 ] 


3 1 00 

as ~~ v °oo3- ~r v ooof 1 ~ x~r S «(« - l)/» 0 -P«-i(^) 

1 d/ ^0 ^0 n = 3 


/d/i 


— |^ooo(l -V) + i I «(« - \)I„ 0 CM 

l 1 n = 3 


The condition dTo/dfl = 0 at 0 = 0 has been used in the above solution. Since if/ = 0 on the bubble 
surface, g(if/) is merely a constant. For the large Re case, it can be established that 


A 0 = 3, /(//) = ( 1+H) 2I \ ^ = 1+ ^"^ ln48 


[26] 


The value of g is obtained by connecting the solution for T 0 from [23] to that given by [20] along 
the front stagnation streamline. The solution for T 0 given above satisfies neither [9] nor [10]. The 
result is good only as r-*\ and cannot be expected to satisfy the condition as r-*co. The adiabatic 
condition at the bubble surface cannot be satisfied as the order of the governing equation has been 
reduced. In fact, the outer field exhibits a logarithmic singularity as the bubble surface is 
approached. In order to satisfy the boundary condition at the bubble surface, we shall develop an 
inner expansion near r = 1. dT 0 /89 does not have a singularity near 9 = 0 (in fact it is zero). Thus 
there is no tangential boundary layer near the forward stagnation point. 

It is seen that f(p) has a singularity at p = - 1; hence a boundary layer (really a thermal wake) 
must exist in the angular coordinate near 0=n. The thermal wake is analyzed in appendix B for 
the large Re flow field. It is shown there that the boundary layer thickness is proportional to e. 
However, the singularity is sufficiently weak that integrals needed for evaluating the bubble 
migration velocity can be obtained without analyzing the thermal wake; from [4] it can be shown 
that the thermal wake contributes to the bubble migration velocity only at 0(e 2 ). 

3.2. The inner temperature field 

In order to satisfy the adiabatic boundary condition at the bubble surface, we need to include 
the neglected conduction term in the radial direction. This is accomplished by suitably stretching 
the radial distance from the bubble surface in an inner (boundary layer) solution. One might 
wonder whether 9 -conduction is also necessary near the forward stagnation point to relieve the 
singularity in the outer solution. The answer is negative. It will be seen from the solution in this 
section that radial conduction alone is sufficient. The physical interpretation is that fluid particles 
on the bubble surface do not infer their temperatures via 9 -conduction from the forward stagnation 
point. Rather, they infer their temperatures via radial conduction across the boundary layer from 
streamlines adjoining the stagnation streamline. The fluid velocity is non-zero everywhere on these 
streamlines and from the outer solution, the temperature along them is bounded. 

We define inner variables (x, p) via x = (r - l)/e and write the inner temperature field t(x, p) 


as 


t(x,p) = t 0 (x,p) + o( 1) 


The energy equation at leading order is 


(du 0 A dt Q .. 2 , 1/2 n ,dt 0 _d% 

v ^ + [lfr ( 1 ’ M y X dx M) Vo(lM) Sp dx 2 

At x = 0 = 0; at 0=0 ^ = 0 

dx o9 


[27] 

[28] 
[29] 


and as x-+oo, the inner temperature field, must match the outer field. A transformation is made 
from (x, p) coordinates to (rj, <£) coordinates where ^ 

[30] 


rj =(1 ~/i 2 ) 1/2 f;o(l, jti)x 




2\ 1 * n(n - 1) 


3 3y 2 (2w — 1) 


LXC n + fp)-C n _ x (p)) 


[31] 
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(£0/) is such that £'(n) = -v 0 (\, n)(\ - n 2 ) 112 , £(1) = 0). rj is related to the streamfunction. This 
transformation reduces the energy equation given above to a diffusion equation with a sink term. 
It is straightforward to find a solution for t 0 as 

f„ = C, + In (f(p)) - i In ( - o + T in ( + 1 F(0 [32] 


where 


C = 


n 

Ul’ 


m= 


D(w) dw 



dw 


[33] 


C, and if are arbitrary constants that will be determined by matching with the outer solution. D(w) 
is Dawson’s function (Abramowitz & Stegun 1968). 

3.2. 1. Matching. The result for the outer field T 0 from [23] written in inner variables, truncated 
at 0(1) in e and rewritten in outer variables is unchanged. The inner result for t 0 from [32] written 
in outer variables and truncated at 0(1) is given below (we use the result that for large C, 
D(C)->1/2C, and hence F(C)->|lnO 


t 0 =C x + In ( fin )) + 




ln(-0 + ^ln 



[34] 


Thus matching requires that 


K = A 0 [35] 

Ci = In e + C 2 [36] 

Even though the inner expansion was initially assumed to be of the form given in [27], the lack 
of a In e term in the outer result forces the choice of C x in [36] to include (l/A 0 ) In e. C 2 is an 0(1) 
constant, the value of which cannot be established unless the outer solution near r = 1 is known 
in more detail, including all the constants. It is straightforward to establish that at 0(ln e), the only 
contribution to the inner field is (l/^4 0 )lne. 

The appearance of In e at leading order in the inner solution requires discussion. A uniformly 
valid expansion can be constructed by additive composition from the inner and outer temperature 
fields. Since the outer field we have given ([23]) is good only near r = 1, the composite expansion 
in this region is merely the inner field. The temperature on the surface of the bubble becomes 
singular (negative infinity) in the limit of e-»0. Of course, in physical systems there is always some 
non-zero conductivity. Therefore the surface temperatures will remain finite. 

The reason for the relatively cold conditions prevalent at the surface of the bubble for small e 
can be understood by examining the behavior of a thin bundle of fluid elements surrounding the 
front stagnation streamline. As the fluid elements approach the bubble from far away, we infer from 
[17] that in the outer region they lose heat and get colder owing to the presence of the sink (the 
i^o term). The reduction in temperature is proportional to the transit time of the fluid elements 
since the sink is constant. This is especially pronounced near the front stagnation point where the 
velocity approaches zero and the transit time diverges logarithmically. When the elements in the 
bundle reach the inner region, they communicate the temperature information by conduction to 
the elements on the bubble surface. Consequently the fluid elements on the bubble surface become 
cold as well. 

The correctness of the above results is qualitatively supported by isotherms sketched from the 
numerical solutions of the problem (Shankar & Subramanian 1988; Balasubramaniam & Lavery 
1989). It is seen that the scaled temperatures on the bubble surface decrease as the Marangoni 
number is increased. The increase is gradual as one might anticipate because of the In e dependence. 
However a precise lne behavior cannot be established from the numerical results as the 
computations have not been performed beyond Ma = 1000. 
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3.3. I„ calculations 

The unknown Fourier coefficients at leading order I„ o , n ^2 are calculated using [4]. It can be 
shown that 


dF \_dto p\ _ ^ooO I 

dn (0, ^~ Z' + 2 A 0 Z 

where £(//) is given by [31]. Thus the Fourier coefficients are given by 


^"0 


CMI—+ 


C 2 


d(i n ^ 2 


[37] 

[38] 


This constitutes a non-linear system of algebraic equations for the unknown Fourier coefficients 
and must be solved numerically, after truncation of the infinite series at a large but finite value 
of n (say n = N). 

A very similar analysis was used by Crespo & Jimenez- Fernandez (1991, 1992) to calculate the 
migration velocity in the limits of Re-+oo and Re->0, for large Marangoni numbers. Our solution 
for large Re is presented below and agrees with their solution. The results for Re->0 are provided 
in section 3.3.2. 

3.3.1. Results for the large Reynolds number flow field. For Re->oo, the inner temperature field 
may be written as 


1 71 1 .2 

t 0 (x,6) = - In £ + 1 + j-j= - - In (27i^ 0 ) + - In 

( = v -fO-nm + n) 


1 + COS Q\ 1 g 2 - 

~toT +«<“«+ jfXC) 


C = - V - rx> l x sin 2 9 

wz 

F( 6= I V(y)dy 
Jo 

is calculated from [5] and [37] as 

Voo0 = 

which computes to 

v coo = t — x In 3 = 0. 1 96007 


1 rv ^ + iv in 30 d0 


4 


o\ Z' 6Z 

l l 


3 8 


[39] 

[40] 

[41] 

[42] 


[43] 


[44] 


Note that [43] is the condition from energy dissipation arguments ([7]). Crespo & Jimenez- 
Fernandez obtain the same equation when they require the mass flux in the 9 direction in the 
momentum boundary layer to be bounded at 9 — n. 

The above result for v w is valid for Re» 1 and Ma» 1. Even though potential flow velocities have 
been used within the thermal boundary layer ([28]) for large Re, there is no restriction that the 
Prandtl number be small. This is because the momentum boundary layer is a layer where the shear 
stress (and not the velocity itself) undergoes a sharp change. As shown in appendix A, the 
correction to the potential flow velocity in the momentum boundary layer is 0(Re -l/2 ). The leading 
order flow field everywhere (including the thermal boundary layer) is the potential flow field 
regardless of the Prandtl number. 

3.3.2. Results for the small Reynolds number flow field. Since the velocity field in the limit Re->0 
is expressed in the form of an infinite series, the Fourier coefficients 7„ 0 cannot be calculated 
analytically and are determined numerically. The infinite series is truncated at a large value of 
n = N. Consistent results for v ot0 have been obtained in the range 10 ^ N < 75. The computations 
were carried out using successive approximations (i.e. fixed point iteration). Basically, a set of 7„ 0 
is guessed and used in the right hand side of [38]. This then generates a new set of I„ 0 and the process 
is repeated till convergence is achieved to a specified tolerance. Convergence is checked for I 2q and 
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the ratio of the sum of the absolute values of the deviations in I„ o in consecutive iterations to the 
sum of the current absolute values of I % . Underrelaxation was used in obtaining the new set of 
/ The required integrals were calculated by Romberg integration (Burden et al. 1981). Typically 
seventeen Romberg steps were necessary for accuracy to four significant digits for the I„ Q . 

It turns out that the computations are quite cumbersome. In fact, we discovered that direct use 
of [24] for the constant A 0 (A 0 is recalculated once a new set of I„ 0 is obtained) is the cause of severe 
difficulties in the calculations. First, odd values of N do not yield a converged set of 7„ 0 by a fixed 
point iteration method; second, for even values of N between 10 and 150, the numerically 
determined values of both A 0 and v ^ diverge with N (they increase monotonically). 

An issue to investigate in connection with the convergence difficulties mentioned above is that 
the temperature gradient on the bubble surface given by [37] is singular at the rear pole ( 9 = n). 
The left-hand side of the energy equation ([8]) has a sink term which, along with the convective 
transport terms, must be balanced by conduction near stagnation points where the velocities go 
to zero. Conduction in the thin radial boundary layer accommodates the sink term at the forward 
stagnation point ( 9 = 0). At 6 - n, the radial boundary layer is infinitely thick; consequently, there 
must exist a boundary layer in the angular coordinate near 9 = n in order to balance the sink term. 
Since the conduction term in the angular direction drops out at leading order in [28] for the inner 
temperature field, the surface temperature gradient given by [37] must be singular at 9 = %. 
However, in the integrals for I„ 0 , this gradient is weighted by CM; it can be shown that 
C n (fi)dT/dn on the bubble surface is bounded everywhere. Thus, in spite of the singularity in dT/dfi 
at 9 =n, it should be possible to determine a set of 7„ 0 . This is precisely what happens for the large 
Re flow field. 

As alluded to earlier, the source of the numerical difficulties is in the expression for A 0 ([24]). 
The series £" =3 h(« - 1)7„ 0 , and more generally the series £“ =3 «(« - 1 )I„ 0 P n -M, appear to have 
convergence problems. The constant A 0 and another constant related to it, viz. 
B 0 = -3y w0 - 5 £"= 3 (- l)" _1 n(« - 1 )I„ are proportional to the viscous normal stress (represented 
by the series given above) at 9 = 0, n respectively. We find that the series for B 0 is divergent. Thus 
the singularity of the normal stress on the bubble surface at 9 = u is controlling the convergence 
of the series in the A 0 expression. 

It is found that if A 0 is taken to be fixed during a fixed point iteration for I„ 0 (i.e. A 0 is not 
recalculated after a new set of I„ 0 is obtained), then the convergence problem for odd N mentioned 
above disappears. Furthermore, typical values of 7„ o thus obtained (non-converged, i.e. not 
summing up to the guessed value of A 0 ) are alternating in sign. Thus the series in [24] for A 0 is 
one of alternating sign and Euler’s transformation (Meksyn 1961) can be used to improve the 
convergence of the series. The Euler summation procedure is used with Van Wijngaarden’s 
algorithm (Press et al. 1986) for calculating A 0 (A 0 recalculated for every new set of 7„ 0 ) together 
with underrelaxation in the fixed point iteration procedure (previous iterates weighted 60-90%). 
This yields a convergent set of I„ 0 (with respect to N) for N in the range 10-75, both for odd and 
even N. The calculated values of 7„ o are presented in table 1 . The final results are 

^0 = 0.1538 [45] 

A 0 = 2.406 [46] 

It is worth making a comment regarding the integrand in the r.h.s. of [38] at the end point 9 = n. 
At first sight, it appears that dT/dn = - V<X)Q /(4B 0 (1 + n)) at this end point. However, B 0 is 
unbounded as noted earlier. As will be clear from a discussion of the nature of the singularity in 
dT/dn to be presented a little later, the correct form of the singularity is dT/dpi 1 +^)~ 3/4 at 
9 = n. Consequently the integrand in [38] is zero at 9 = n. We have used this knowledge in 
calculating the values of the constants 7„ 0 reported in table 1. 

Streamlines in a laboratory reference frame are displayed in figure 1 . The topology of the flow 
field is virtually unchanged from that shown by Shankar & Subramanian (1988) for Ma = 200. The 
saddle point behind the bubble has moved closer to the bubble, when compared to their plot. 
Figure 1 clearly shows the transition from a potential dipole dominated flow field ahead of the 
bubble driven by the P, (Legendre) mode of the surface temperature distribution to that driven 
by the P 2 mode behind the bubble (see discussion in Shankar & Subramanian). 
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Table 1. The calculated values of the Fourier coefficients for the small Re 

flow field 


n 

/„ 

n 

4 

2 

-0.30761 

14 

-0.000747 

3 

0.06087 

15 

0.000624 

4 

-0.02557 

16 

-0.000528 

5 

0.01237 

17 

0.000451 

6 

-0.00750 

18 

-0.000389 

7 

0.00477 

19 

0.000339 

8 

-0.00334 

20 

-0.000297 

9 

0.00241 

21 

0.000262 

10 

-0.00182 

22 

-0.000233 

11 

0.00141 

23 

0.000208 

12 

-0.00112 

24 

-0.000187 

13 

0.000906 

25 

0.000169 


The temperature drop over the bubble surface from the front pole to the rear pole may be shown 
to be 

CO 

A7- = - £ (4/> - 1 )I 2p [47] 

P= > 

for the low Re flow field. This relation is obtained by integrating the tangential stress balance at 
the bubble surface ( d/dr)(v/r ) = dT/86 at r = 1. From the form of the singularity mentioned above, 
it is seen that (8T/86)( 1, 6) at 6 = n is integrable. Thus the series in [47] is expected to converge. 
Since the infinite series in [1] to [3] have been truncated at n = N for numerical evaluation of /„, 
the series in [47] is summed to the same level of truncation. Numerically determined values of AT 
for various N obtained in this manner are shown in figure 2. It is interesting that AT plotted versus 
1 is fitted well by a straight line and extrapolates to a value very close to n/2 as iV->oo. 



Figure 1. Streamlines (equally spaced in !F) in a laboratory reference frame for the case Re = 0, Ma-+oo. 
The values of the Fourier coefficients used to calculate the streamlines are given in table 1 . 
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Figure 2. Scaled temperature difference, AT, over the bubble surface vs l/^/iv. N is the truncation level 
of the infinite series in [1]. The symbols (□) represent the actual values from numerical calculations. The 
straight line is a linear regression line. The intercept which represents AT for N-+ oo is 1.5715 + 0.00617 

which is very close to n/2. 


We discuss below the nature of the temperature singularity near the rear stagnation point. It must 
be stated at the outset that this singularity is unimportant in determining I n0 ; the use of the Euler 
summation procedure has circumvented the convergence difficulties and enables the calculation of 
7„ 0 . This thermal singularity is produced because 9 -conduction has been neglected in the energy 
equation near 6 = n. The scalings (i) u • VT ~ ~ 1 near the rear stagnation point and (ii) shear 

stress ~VT(it turns out that for low Re the components of the stress tensor, in particular dv/dO, 
also scale as |Vr|) can be shown to yield v ~ [VZ’j “ 1 ~ (n — 0) 1/2 ~ (1 + /r) ,/4 . 

It is also possible to infer the nature of the singularity in dT/dji from the calculated set I„ Q , In 
order to do this, we construct an analog of the Domb-Sykes plot used in the analysis of singularities 
of functions represented by power series (see Van Dyke 1975). Consider first the model function 
1/(1 + n) a for 0 <a < 1. This can be expanded in a series of Legendre polynomials as follows 

a irf™ [481 


where 


H n = {— 1)"(2« + \)2~ a 


r(l-fl) r(n+a) 
r(a) r(n + 2-a ) 


[49] 


A plot of versus \/n yields a straight line with the slope (1 — 2d) for large n. The scaled 

viscous stress on the bubble surface can be written as 


S(fi)= \ Z «(» “ W no P n -M [50] 

By constructing a plot similar to that for the model function, one might guess that the singularity 
is the same as that of the model function. Such a plot is shown in figure 3 in the region where the 
ratios of the Fourier coefficients appear to be settling down after displaying oscillatory behavior 
for small values of n. The slope is 0.6 which gives a value a = 1/5, in contrast to the exact value 
a = 1/4. We believe that a greater precision is needed in I„ 0 to extract the correct value of a from 
the numerical calculations. 
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1/n 


Figure 3. Plot of the ratios of successive Fourier coefficients vs 1/n for the infinite series in the expression 

for the normal stress in [50]. 

For the large Re flow field, the tangential stress balance cannot be satisfied without including 
the momentum boundary layer. Therefore a relation for A T analogous to [47] cannot be obtained. 
Also, the singularity in dt/d6( 1, 9) at 9 = n from [39] is not integrable. Hence the boundary layer 
in the angular direction near 9 = n has to be analyzed carefully to determine AT over the bubble 
surface. This issue is addressed in appendix B. 


4. HIGHER ORDER THERMAL BOUNDARY LAYER ANALYSIS 

The leading order outer and inner temperature fields (T 0 and * 0 ) were analyzed in section 3. The 
higher order corrections in e for the large Re flow field are calculated in this section. The chief 
assumption is that while Ma and Re are both large, the Prandtl number is small. That is, potential 
flow velocities are used in the convective terms in the energy equation and the viscous boundary 
layer contribution to energy dissipation is ignored. The consequence is that effects of order 
0(Re -l/2 ) are neglected. It is possible to obtain the proper 0(e) corrections to the temperature field 
since inclusion of the momentum boundary layer effects will yield corrections of O(ePr) which will 
be negligible for Pr« 1 . Hence the thermal corrections to the migration velocity are more important 
than corrections from the inclusion of the momentum boundary layer effects. As will be evident 
later, the analysis provides corrections to the bubble migration velocity at orders 0(e lne) and 
0(e). 

4.1. Outer temperature field 

Let the outer temperature field be 

T=T 0 + eT l + o(e). [51] 

Here T 0 is the leading order outer temperature field given by [20]. Examining the governing 
equations for T, it can be verified that T, satisfies homogeneous equations and boundary conditions 
and is thus identically zero. Thus non-trivial corrections to T 0 occur only at o(e). In section 3, T 0 
was obtained via Taylor series expansions of the velocity fields in (r — 1). In this section we 
construct solutions to T 0 formally that are good to 0(e) when rewritten in inner variables. This 
is necessary in order to extend the inner solution to that order. 
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We note that T 0 may be expressed alternatively as 


T 0 (r,9) = GW o)- 


2 r A 


2r 3 + 1 


d 0 
sin 0 


[52] 


where the integrand in the indefinite integral must be expressed as a function of i// 0 and 0. Since 
ijj 0 is a parameter in [20] and [52], one can expand T 0 about a given value of i J/q, in particular i// 0 = 0. 
We use [20] to expand T 0 in the vicinity of 0 = 0 and arbitrary r and [52] to expand T 0 in the vicinity 
of r = 1, 0 ^ 0 ^ n. These two expressions for T 0 are then “patched” in the neighbourhood of r = 1 
and 0=0 so that they are representations of the same function. This process yields 
G(il/ 0 ) = a 0 + b 0 \n 11/0 + ^11/0 + b^ 0 \nil / 0 and the values of the constants are determined by the 
patching procedure mentioned above. After a lengthy but straightforward calculation we may 
determine T 0 near r = 1 (up to O(ij/ 0 )) via [52] to be 


«'■ 9) “ (‘ + ijf" 432 ) ~Ts (fi + ln432 )( r2 ~ 9 ^ e + 5 ,n ( r2 ^ 

+ ? In (1 + cos 0) + ^ r 2 - ^ cos 0 + ~ sin 2 9^r 2 - ^ In ^ 


+? sin 2 6^r 2 - ^ In (1 + cos 0) 
We may now write T in inner variables as 


[53] 


T(x, 0) = ^ln^ + ^l+-^/|-^ln48 + ^lnx + ^ln(l + cos 0)^ + e In e^-x sin 2 0^ 


'2 1 2 

-\-e{-x cos0 +-x sin 2 0 \nx +-x sin 2 0 in (1 + cos0) 
'3 3 3 


1/ n 


6^3 


+ In 48 lx sin 2 0 


[54] 


Thus the outer temperature field has correction terms at orders e In e and e when written in inner 
variables. This demonstrates that the next corrections to the inner field occur at 0(e In e) and O(e). 
These are calculated below. 


4.2. Inner temperature field 
Let 


Voo = v x0 + e In e v^ n + ev xl + o(e) 


[55] 


t = i In e + t 0 + e In e t n + et x + o(e) [56] 

t 0 has been obtained in section 3 ([39]). In what follows, we will obtain solutions for t tl and t\ and 
also calculate v^n and v^i . We use the match with [54] to establish boundary conditions on t^ and 
as x->oo. 

4.2.1. Inner field at 0(e lne). From [8], the equation for t n may be written as 


^coo 


. dt\\ 3 . _ dt [ i 
- 3x cos 0 ~ + - sin 0 — 
ox 2 cv 


d\ 

dx 2 


Vooii 

v<x)Q dx 2 


[57] 


The boundary conditions are 


d fi-Q atJc-0 
dx 


[ 58 ] 
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1 • 20 

xsnr0 asx->oo 
3 

[59] 

— = 0 at 0 = 0 
50 

[60] 


Equation [60] and a similar condition on t x follow from symmetry. The solution is 

’ = 'ft D ' ( ° + ^~ 2Vl ( ierfC C " e i) 1611 

where ip = (§Ko* sin 2 0, £ = ( 5 K 0 O - cos 0) 2 ( 2 + cos 6) and C = $ should really be 
^ 0 , but we omit the subscript in this section for convenience. 

The contribution to the bubble velocity at 0(e lne) may be calculated from [7], [55] and [56] 
as 


^co/l 


1 

4 


'n 

Jo 



sin 2 0 d6 


This yields 


^oo/l 


sin 5 6 


72jc» m0 Jo (1 -cos0)(2 + cos0) ,/2 


d9 


= - -0.1369 

V 72^0 \ 7 21 / 


4.22. Turner field at 0{e). The equation for can be obtained from [8] to be 

J/, , 3 . 5f,\ 5 2 t, 5 / 5? 0 \ ^,5% 

2x — 3x cos 0— — h - sin 0 — = -r— j + 2 — I x— 1 -j-j 

dx 2 50 / 5x 2 5x V 5 a: 7 v m0 ox 


with boundary conditions 


dt x 

dx 


= 0 at x = 0 


^ = 0 at 9 = 0 

50 


[62] 

[63] 

[64] 

[65] 


and ti matches T 0 ([53]) far away. The solution can be found as the superposition of a particular 
solution and three homogeneous solutions t m , t xh2 , that are constructed below. 

• The particular solution is 


- H cos ® + sin2 6 in (H^)) ■ ■ 2cm)) 

+ w(6)(-2C + (4£ 2 - 2)D(Q) [66] 

where 

HO) = jfri ^ cos 0 - i cos 2 0 + ^ cos 4 0 — ^ [67] 

• At \]/=0 dt m /dij/ = 0; as ^->oo J- 71/(2^) -Jin 48V; at 0=0 

dt m /d9= 0. The solution is \ V*W J 


t\h\ ~ 


9v a 


In 


3v 


ooO 


~—p — 1 In 48 

2^3 2 


)(^-2Vj(cerfcC ' *-«>)) 


[ 68 ] 
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At ij/ =0 t m = 0; as 1/7-00 (2/9(i; oo0 ))^ lnft/r/2). The solution is 

't 

9fooo V 2 




where 


(0/y + (l + 2C 2 )L' = 2£ 


T(0 = 


' C P -Dip) 


dp. 


[69] 

[70] 

[71] 


4/ w 

At ij/ = 0 dt m /di{/ = -(4/(917*0)) (cos 0/sin 2 0 + In ((1 + cos 0)/sin 0)) + 2 w(8)/y/t - (d/, M / 
d\f/ )^, = 0 = AT(<^) + A:/V^; as ifr -oo r 1 /l 3 -0; At 0 = 0 3^ 3 /S0 = 0. A: is chosen such that tf(0 

is not singular at 0 = 0 and is given by 


27 


8 u r 


The fundamental solution that satisfies the condition dP/d\j/ = 1 at if/ = 0, P—0 as 1/7 — 00 and 
5P/00 = 0 at 0 = 0 is 


/*(*, o = Vf (c erfc t - e ' !2 j 

/ 1W can be obtained from Duhamel’s theorem (Myers 1987) to be 

r- dP 

t m (&, O = —ky/n erfc C + I K(z) -(&,£- t) dz 

1 *■ 


[72] 


exp - 


4(£ -t) 


dt 


[73] 


The contribution to the bubble velocity at this order, 0 *,, is determined via [7], [55] and 
[56] to be i 7 *,= -iJS(0^/00)(l,0)sin 2 0d0. It can be shown that at r = 1, 
(dtJdO) = (dt m fd6) + (dt m /d6). From [73], it can be determined that / 1A3 (1, 0) has a logarithmic 
singularity near 0 =jc. After integrating the expression for 0 *, by parts we finally obtain u*, to 
be 


1 


V„, = 


' n r° /\ 4 - cos a\ + cos a ^ 3 



0 jo 


(3 + cos a) 


sin a I sin 2 a 4 (1 — cos a)(2 + cos a) 2 


3y/n 

1 9^0 \ sin 0 cos 0 sin 3 a , 

+-lm + p — 7 — da d 0 

4 4^/r / V / <^(6)-T(a) 


— 0.1369( In. „ 

\ V 3 y oo 0 


-V-^= + ln48 

2 VV3 


[74] 


where t(a) = ( 17 * 0 / 2 ) (1 - cos a) 2 (2 + cos a). The double integral is evaluated numerically. The 
result is 

ti*! = 0.5311 +0.1267 = 0.6578 [75] 

Thus for the case Re— 00 , the result for the bubble migration velocity up to 0(e) is 


v * = ( - - — ) - 0.1369c In £ + 0.6578c. 

CO l 3 g I 


5. DISCUSSION 


[76] 


Crespo & Jimenez-Fernandez (1992) report a leading order result v oo0 = 0.223 for the case of 
Re— 0, They mention using N = 20 for their calculations and report good convergence of the result. 
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While details of their numerical procedure to calculate the Fourier coefficients are not given, we 
suggest that the value reported by them is not correct. To establish this, we performed calculations 
for N - 20 with A 0 being calculated at every iteration without the benefit of the Euler summation 
and obtain v^o = 0.2289, A 0 = 4.452. y TO0 is close to the value obtained by Crespo & Jimenez- 
Fernandez (note: Crespo & Jimenez-Fernandez do not report the value of A 0 and therefore this 
cannot be checked). In view of our comments in section 3.3.2, we do not believe this to be the 
converged value of v^. 

To establish the internal consistency of our results, we checked the equality of the left and right 
sides of the following equation obtained by integrating the energy equation ([8]), where [4] for /„ 
has been inserted 

1 co Ip q2 T 

^(l -30 = 4 “ 1)/ ' + M^J _ 1 

i co i p r 1 

1,11 

This equation may also be obtained alternatively from the viscous dissipation arguments, 
mentioned in section 2, applied to the low Reynolds number velocity field. Our numerical results 
satisfy the above equation. 

A key result from the present leading order analysis is that , the scaled bubble velocity, reaches 
an asymptotic limit for large Marangoni numbers. The limit for Re->0 is smaller than that for 
Re->oo by 27.4%. This supports the behavior determined numerically by Balasubramaniam & 
Lavery (1989), who found the bubble velocity for Ma = 1000 to increase with Re, though no 
asymptotic behavior was evident. 

Shankar & Subramanian (1988) report that their scaled bubble velocity from numerical 
calculations under the creeping flow assumption decreases logarithmically with the Marangoni 
number. The extrapolated value from their curve fit to the numerical data is zero as the Marangoni 
number goes to infinity. We conclude that the behavior they predict is only true locally in the range 
of Ma in which the results were curve-fitted (75 ^ Ma ^ 200). 

The higher order analysis performed in section 4 shows that, at least for the large Re flow field, 
logarithmic dependence of the form e In e occurs in the bubble velocity. Thus the approach to the 
asymptotic velocity is quite slow and the limit is attained only for very large Marangoni numbers. 
This explains why asymptotic behavior was not discerned in previous numerical simulations. 

The solution for the thermal wake for large Re and Ma obtained in appendix B is interesting. 
The chief limitation of the analysis is that the use of potential flow velocities in the energy equation 
leads to a temperature drop over the bubble surface of order 0(ln e) which is physically unrealistic. 
When the momentum wake analyzed in appendix A is taken into account in the thermal wake 
region, we speculate that this temperature drop will be of order 0(1). 


6. CONCLUDING REMARKS 

It has been demonstrated here that the scaled thermocapillary migration velocity approaches a 
non-zero asymptotic limit as the Marangoni number Ma~>oo. This is true in the two limits Re->0 
and Re ->oo. Given the nature of the results, it would be reasonable to conjecture that the statement 
holds true for any Re [see Crespo & Jimenez-Fernandez (1992) who correctly comment that a 
similar analysis of the thermal boundary layer could analytically provide a boundary condition for 
the flow field for arbitrary Re]. We have also shown that the summations involved in the low 
Reynolds number problem must be handled carefully in order to arrive at the correct result for 
the migration velocity. 

The analysis has also been extended to calculate the higher order corrections to the migration 
velocity for the case Re->oo. Our principal result is that the next correction occurs at 0(e lne) 
and not at 0(e) as one might guess. While a similar extension is possible for Re->0, the calculations 
for the Fourier coefficients appear formidable and we do not attempt them in this work. 

The thermal wake has been analyzed for the case Re->oo. In principle, the results can similarly 
be extended for Re->0 as well. The thermal wake does not influence the bubble velocity up to 0(e). 
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APPENDIX A 


The Momentum Boundary Layer for Large Reynolds and Marangoni Numbers 

The temperature field calculated in section 3 for the large Reynolds number flow field assumes 
that the velocities are given by the potential flow solution. We shall show below that this is a good 
assumption everywhere except near the rear stagnation point. There is a thin boundary layer near 
the bubble surface where the shear stress changes from that in potential flow (outside the boundary 
layer) to that demanded by the tangential balance at the bubble surface. This boundary condition 
is 


at r = 1 


dT 

dr\r) Jd 

From the result in section 3.3.1, the temperature gradient at the bubble surface is 




+ 


sin 3 6 


2(1 — cos 9) 2 ( 2 + cos 9) 


[Al] 


[A2] 


Following Moore (1963), the velocity field is written as the sum of the potential flow field and a 
boundary layer correction field. The momentum boundary layer problem is the same as in Moore’s 
analysis, except for a change in the shear stress boundary condition. 


8v Vdv IdV _J_^v 
U dr r 89 r 39 V Re dr 2 


[A3] 


Here v is the 9 -component of the velocity correction in the boundary layer; Re = ( — o r )AR 2 /(pv) 
is the Reynolds number; XJ, V are the potential flow velocity components near the bubble surface. 
Using suitable Taylor series expansions for U and V, the boundary value problem can be written 
as 


, dv 3 . „ dv 3 

^ooo( — 3(r-l)cos0— + -sm9~+-vcos9 

v — 0 at 9 — 0 and as r -> oo 

TP = 3^ sin 9 + T' s (9) at r = 1 
dr 


\__d_h 
Re dr 2 


[A4] 
[A 5] 
[A6] 


The dependent and independent variables are transformed as follows: 

p = ^/Re v sin 9, x =|\/Re^<»osin 2 0(r — 1), £ = ^(1 — cos 9) 2 (2 + cos 9) [A7j 


dp_d 2 p 

dd, ~ d X 2 


p = 0 at £ = 0 and as x -+ °o 


dp 

h 


= <?«)« 2 - 


4 

91^ sin 2 9 


+ 


sin 2 9 


This yields 


at rj = 0 


[A8] 

[A9] 

[A10] 
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It may be noticed that £ is a time-like variable. The above “initial” value problem can be solved 
using Duhamel’s theorem (Myers 1987). The fundamental solution, P(£, /), which satisfies the 
boundary condition dP/dx = 1 at x = 0, is the same as Moore’s solution and is 


P(£, x ) = Ce rfcC-4= e_C2 

71 


C = 


2 Jl 


Thus by Duhamel’s theorem 

p(£,x) = 


i dP If- 

G(t)— (£ -T,x)dT = 7= 

Jtij o 




Gix) 


exp - 




dr 


[All] 


[A 12] 


The expression for G(g) ([A 10]) is complicated and the Duhamel integral has to be evaluated 
numerically. Integrating />(£, x) with respect to x> the result for the volumetric flow rate in the 
boundary layer given by Crespo & Jimenez-Fernandez (1991) can be recovered. 

p (£ x ) given by [A12] is expected to be of order 0(1); consequently v(r, d ) in the momentum 
boundary layer scales as Re 1/2 . Thus for large Re, the potential flow is only slightly perturbed in 
the momentum boundary layer [0(1) changes occur only in the velocity gradient, i.e. stress]. It is 
therefore justified to include only the potential flow portion of the velocity field in the energy 
equation to calculate the temperature field. This is questionable near the rear stagnation point, 
where v from [A 12] is expected to be singular. 


APPENDIX B 


The Thermal Wake for Large Reynolds and Marangoni Numbers 

In section 3, it is seen that both the outer and inner temperature fields are singular at 6 = n. This 
happens for both the large and small Reynolds number flow fields. We shall attempt to remove 
this singularity by analyzing the tangential thermal boundary layer near 0 = it, for the case Re->-co 
using the potential flow velocities given in section 2. 

The potential flow velocities given by [6] are not strictly valid in the thermal wake region. This 
is because the momentum boundary layer analyzed in appendix A will separate and form a flow 
wake in the vicinity of the rear stagnation point. Though the vorticity vanishes on the rear 
stagnation line 6 = n and the flow field remains irrotational near it, the corrections to the velocities 
given by [6] will not be negligible in the wake. The flow wake will have a significant impact on 
the thermal wake. Balasubramaniam (1995) has compared the temperature profile in the thermal 
wake via numerical computations with and without the use of the potential flow velocities from 
[6] (numerically determined velocities are used in the latter case). He concluded that the two thermal 
wakes are drastically different away from the bubble. Nevertheless, for the purposes of this idealized 
analysis, we shall assume that the velocities are given by [6]. 


Outer temperature field 

The temperature field in the absence of conduction must first be determined. This is done in a 
manner similar to the r -outer analysis performed in section 3.1 (we use the notation r -outer and 
“r -inner” to refer to the outer and inner solutions in the radial coordinate). The outer equation 
near 6 = n is 


1 \6T 4> 


1+ 1-A 1+^ ^=0 


7 dr 


1 \dT 


2r 2 J d(p 


[Bl] 


where 0 = n - 0 and 0 «1. As r-*oo, WT must attain the imposed uniform field and near the 
surface of the bubble, T must agree with the result in section 3.1 ([23]). The solution is 

fi 2 + r + 1 


= C 3 + -ln0 -r +-ln r 2 -- +-ln 


1\ 1 


Ci = 2 — 


71 


1 


6^3 


<r " l) 2 

In 432 


1 (2 r + 1\ 

+ 75 arctan W 


[B2] 

[B3] 
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The outer field given above and the r -inner temperature field given by [39] are both singular at 
0 = n. It is convenient to construct an r -composite field T c near 9 — n from these two results. The 
composite field serves as the -outer field for the thermal wake. However, it turns out that T c 
constructed from [39] and [B2] is not valid in the vicinity of </> = 0. Substitution of this T c in [8] 
reveals that the remainder is proportional to (r — l)(r 2 — \)D\£)/r A , which is non-zero [i.e. 0(1)] 
at <f> = 0 (i.e. C = 0) and r not equal to 1 and therefore not tolerable. The source of this behavior 
is the F(C) term in [39]; its contribution to T c is such that u • VT C leaves the remainder given above. 
Fortunately, it is trivial to overcome this difficulty by seeking a solution in terms of the variable 
A given below. F{X) has no contribution to the convective terms (A is proportional to the stream 
function near 6=n; £ is proportional to the stream function near r = 1). The composite 
temperature field may be written as 

+ [B4] 

where _____ 

1= (r 2 --)^ [B5] 

V 32 V V 


Inner temperature field 

The equation for the temperature field in the wake including conduction in the tangential 
direction is (the inner variable is y = </>/e) 


VoQO 


1 + ( 1 — 


dt y 
dr r 


J _\dT 
2r 3 J dy _ 


\ _d_/ dt\ 
r z y dy \ dy) 


[B6] 


with boundary conditions 


_ = 0 at y = 0 
dy 


[B7] 


^ = finite at r = 1 [B8] 

dr 


and matching with T c written above [note that F(X) rewritten in inner variables and truncated to 
0(1) becomes F( 0) which is merely a constant]. An important observation with respect to the 
boundary condition in r is that at leading order t must not contain a term of the form ln(r — 1) 
and thus dt/dr must be finite at r = 1. We cannot impose dt/dr = 0 at r = 1 as this precludes the 
occurrence of higher order boundary layers (see section 4). An alternative boundary condition 
would be that the r -inner temperature field in the wake region [obtained by rewriting t in (x,y) 
coordinates and truncating to 0(1)] must satisfy the adiabatic condition dt /dx = 0 at * = 0. 

A similarity solution to t can be obtained to be 


«(r.)-)-C J + ?F(9) + ilng + ln £ +?ln,+}ln(r» — 7 ) + g H 

1 /2r + l\ l p 

_r + _ arctar \ 73 / 3 ' 


r 2 + r + 1' 

0* - 1?" 




r 2 + r +'l 


[B9] 


where E x (x) = J® (e~ p (p) dp is an exponential integral (Abramowitz & Stegun 1968). 

The composite solution in the wake is the same as t with F( 0) replaced by F(X). It is seen that 
while dt /dr — 1 as r-»oo, t-> — r -}-§ In r. This behavior is tolerable at leading order, but points 
to difficulties that might arise at higher orders. 

It can be shown that the temperature difference over the bubble surface from the forward to the 
rear stagnation points scales as — jin e, and is thus unbounded as Ma->co. This agrees with the 
findings of Balasubramaniam (1995) who pointed out that this is a consequence of using the 
potential flow velocities in the energy equation. Indeed for the Re->>0 flow field, AT over the bubble 
surface is finite ([47]), while for large Re, numerical calculations (Balasubramaniam & Lavery 1989) 
do not seem to indicate that ATis unbounded. It remains to be shown for large Re that when proper 
account is taken of the velocity corrections in the momentum wake, AT is indeed bounded. 



